#ifndef MATRIX_H
#define MATRIX_H

#include <QDebug>
#include <QObject>
#include <QTime>

#include <Eigen/Dense>
#include <complex>

#ifdef IMPL_MATRIX_OPERATOR
#undef IMPL_MATRIX_OPERATOR
#endif
#define IMPL_MATRIX_OPERATOR(TYPE)     TYPE& operator +(const TYPE& other)   \
                                       {                                     \
                                            m_data = m_data + other.m_data;  \
                                            return *this;                    \
                                       }                                     \
                                       TYPE& operator -()                    \
                                       {                                     \
                                            m_data = -m_data;                \
                                            return *this;                    \
                                       }                                     \
                                       TYPE& operator =(const TYPE& other)   \
                                       {                                     \
                                            m_data = other.m_data;           \
                                            return *this;                    \
                                       }                                     \
                                       TYPE& operator -(const TYPE& other)   \
                                       {                                     \
                                            m_data = m_data - other.m_data;  \
                                            return *this;                    \
                                       }                                     \
                                       TYPE& operator +=(const TYPE& other)  \
                                       {                                     \
                                            m_data += other.m_data;          \
                                            return *this;                    \
                                       }                                     \
                                       TYPE& operator -=(const TYPE& other)  \
                                       {                                     \
                                            m_data -= other.m_data;          \
                                            return *this;                    \
                                       }                                     \
                                       TYPE& operator *(const Scalar& scalar)\
                                       {                                     \
                                            m_data = m_data * scalar;        \
                                            return *this;                    \
                                       }                                     \
                                       TYPE& operator *(const TYPE& other)   \
                                       {                                     \
                                            m_data = m_data * other.m_data;  \
                                            return *this;                    \
                                       }                                     \
                                       TYPE& operator /(const Scalar& scalar)\
                                       {                                     \
                                            m_data = m_data / scalar;        \
                                            return *this;                    \
                                       }                                     \
                                       TYPE& operator *=(const Scalar& scalar)\
                                       {                                     \
                                            m_data *= scalar;                \
                                            return *this;                    \
                                       }                                     \
                                       TYPE& operator /=(const Scalar& scalar)\
                                       {                                     \
                                            m_data /= scalar;                \
                                            return *this;                    \
                                       }

typedef std::complex<float> ComplexFloatType;
typedef std::complex<double> ComplexDoubleType;

Q_DECLARE_METATYPE(ComplexFloatType)
Q_DECLARE_METATYPE(ComplexDoubleType)

// TODO: add reference count
class MatrixBase : public QObject
{
    Q_OBJECT
public:
    enum Type { Int, Float, Double, ComplexFloat, ComplexDouble };

    MatrixBase(int rows, int cols) { Q_UNUSED(rows) Q_UNUSED(cols) }
    explicit MatrixBase() {}
    virtual ~MatrixBase() { qDebug() << "now releasing matrix data ..."; }

public:
    virtual Type type() const = 0;
    virtual int columns() const = 0;
    virtual int rows() const = 0;
    virtual void random() = 0;

    virtual void transpose() = 0;
//    virtual void conjugate() = 0;
    virtual void adjoint() = 0;

    virtual QString toJSON() const;
    void serialize(QDataStream& out) const;
    void unserialize(QDataStream &in);

signals:
    void valueChanged(int row, int col, const QVariant& oldValue, const QVariant& newValue);

protected:
    virtual QString toJSONHelper() const = 0;
    virtual void serializeHelper(QDataStream& out) const = 0;
    virtual void unserializeHelper(QDataStream& in) = 0;
};

template <typename T>
static inline bool isEqual(const T& lhs, const T& rhs)
{
    return lhs == rhs;
}

template <>
static inline bool isEqual<float>(const float &lhs, const float &rhs)
{
    return qFuzzyCompare(lhs, rhs);
}

template <>
static inline bool isEqual<double>(const double &lhs, const double &rhs)
{
    return qFuzzyCompare(lhs, rhs);
}

template <>
static inline bool isEqual<std::complex<float> >(const std::complex<float> &lhs, const std::complex<float> &rhs)
{
    return qFuzzyCompare(lhs.real(), rhs.real()) && qFuzzyCompare(lhs.imag(), rhs.imag());
}

template <>
static inline bool isEqual<std::complex<double> >(const std::complex<double> &lhs, const std::complex<double> &rhs)
{
    return qFuzzyCompare(lhs.real(), rhs.real()) && qFuzzyCompare(lhs.imag(), rhs.imag());
}

template <typename T>
class Matrix : public MatrixBase
{
public:
    typedef typename T::Scalar Scalar;

    Matrix(int rows, int cols) : MatrixBase(rows, cols)
    {
        try {
            m_data = T(rows, cols);
        } catch (std::exception& exp) {
//            throw exp;
            qDebug() << exp.what();
        }
    }

    explicit Matrix(const Scalar* data)
    {
        m_data = T(data);
    }

    Matrix(const Matrix& other)
    {
        m_data = other.m_data;
    }

    int columns() const { return m_data.cols(); }
    int rows() const { return m_data.rows(); }
    void random() { m_data.setRandom(rows(), columns()); }

    void resize(int rows, int cols) { m_data.resize(rows, cols); }

    void setValue(int row, int col, Scalar value)
    {
        Scalar oldV = m_data(row, col);
        if (!isEqual<Scalar>(oldV, value)) {
            m_data(row, col) = value;
            emit valueChanged(row, col, QVariant::fromValue<Scalar>(oldV), QVariant::fromValue<Scalar>(value));
        }
    }
    const Scalar& value(int row, int col) const { return m_data(row, col); }

    void transpose() { m_data.transposeInPlace(); emit valueChanged(-1, -1, QVariant(), QVariant());}
//    void conjugate() { m_data.conjugateInPlace(); }
    void adjoint() { m_data.adjointInPlace(); emit valueChanged(-1, -1, QVariant(), QVariant());}
    Scalar sum() const { return m_data.sum(); }
    Scalar product() const { return m_data.prod(); }
    double mean() const { return m_data.mean(); }
    Scalar minCoeff() const { return m_data.minCoeff(); }
    Scalar maxCoeff() const { return m_data.maxCoeff(); }

    QVector<Scalar> row(int row) const
    {
        int cols = m_data.cols();
        QVector<Scalar> elements(cols);
        for (int i = 0; i < cols; ++i)
            elements[i] = m_data(row, i);
        return elements;
    }

    QVector<Scalar> column(int col) const
    {
        int rows = m_data.rows();
        QVector<Scalar> elements(rows);
        for (int i = 0; i < rows; ++i)
            elements[i] = m_data(i, col);
        return elements;
    }

protected:
    T m_data;
};

class IntMatrix : public Matrix<Eigen::MatrixXi>
{
public:
    IntMatrix(int rows, int cols);

    Type type() const { return Int; }

    IMPL_MATRIX_OPERATOR(IntMatrix)

protected:
    virtual QString toJSONHelper() const;
    virtual void serializeHelper(QDataStream& out) const;
    virtual void unserializeHelper(QDataStream& in);
};

class FloatMatrix : public Matrix<Eigen::MatrixXf>
{
public:
    FloatMatrix(int rows, int cols);

    Type type() const { return Float; }

    IMPL_MATRIX_OPERATOR(FloatMatrix)

protected:
    virtual QString toJSONHelper() const;
    virtual void serializeHelper(QDataStream& out) const;
    virtual void unserializeHelper(QDataStream& in);
};

class DoubleMatrix : public Matrix<Eigen::MatrixXd>
{
public:
    DoubleMatrix(int rows, int cols);

    Type type() const { return Double; }

    IMPL_MATRIX_OPERATOR(DoubleMatrix)

protected:

    virtual QString toJSONHelper() const;
    virtual void serializeHelper(QDataStream& out) const;
    virtual void unserializeHelper(QDataStream& in);
};

class ComplexFloatMatrix : public Matrix<Eigen::MatrixXcf>
{
public:
    ComplexFloatMatrix(int rows, int cols);

    Type type() const { return ComplexFloat; }

    IMPL_MATRIX_OPERATOR(ComplexFloatMatrix)

protected:
    virtual QString toJSONHelper() const;
    virtual void serializeHelper(QDataStream& out) const;
    virtual void unserializeHelper(QDataStream& in);
};

class ComplexDoubleMatrix : public Matrix<Eigen::MatrixXcd>
{
public:
    ComplexDoubleMatrix(int rows, int cols);

    Type type() const { return ComplexDouble; }

    IMPL_MATRIX_OPERATOR(ComplexDoubleMatrix)

protected:
    virtual QString toJSONHelper() const;
    virtual void serializeHelper(QDataStream& out) const;
    virtual void unserializeHelper(QDataStream& in);
};

QDebug operator<<(QDebug, const IntMatrix &);
QDebug operator<<(QDebug, const FloatMatrix &);
QDebug operator<<(QDebug, const DoubleMatrix &);
QDebug operator<<(QDebug, const ComplexFloatMatrix &);
QDebug operator<<(QDebug, const ComplexDoubleMatrix &);


#endif // MATRIX_H
